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ABSTRACT 

We compare N-body simulations of isolated galaxies performed in both frameworks of modified Newtonian dynamics 
(MOND) and Newtonian gravity with dark matter (DM). We have developed a multigrid code able to efficiently solve 
the modified Poisson equation derived from the Lagrangian formalism AQUAL. We take particular care of the boundary 
conditions that are a crucial point in MOND. The 3-dimensional dynamics of initially identical stellar discs is studied 
in both models. In Newtonian gravity the live DM halo is chosen to fit the rotation curve of the MOND galaxy. For 
the same value of the Toomre parameter (Qt), galactic discs in MOND develop a bar instability sooner than in the 
DM model. In a second phase the MOND bars weaken while the DM bars continue to grow by exchanging angular 
momentum with the halo. The bar pattern speed evolves quite differently in the two models: there is no dynamical 
friction on the MOND bars so they keep a constant pattern speed while the DM bars slow down significantly. This 
affects the position of resonance like the corotation and the peanut. The peanut lobes in the DM model move radially 
outward while they keep the same position in MOND. Simulations of (only stellar) galaxies of different types on the 
Hubble sequence lead to a statistical bar frequency that is closer to observations for the MOND than the DM model. 
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1. Introduction 



Galaxies: kinematics and dynamics — Galaxies: spiral — Galaxies: structure 



As has been emphasized in the last years, the concordance 
ACDM cosmological model is very successful in account- 
ing for large-scale structure formation (e.g., Silk 2004), but 
encounters severe problems at galactic scale: in particu- 
lar the highly peaked dark matter (DM) distribution pre- 
dicted by numerical simulations (Navarro et al. 2004) is not 
compatible with most observed rotation curves of galax- 
ies (de Blok 2005); the predicted angular momentum of 
baryons condensed in galaxies is much too low (Steinmetz 
2003), and the number of predicted satellites around a given 
giant galaxy is more than an order of magnitude larger 
than what is observed (Moore et al. 1999). One solution 
to these problems has been searched for in the energetic 
feedback provided either by violent star formation (e.g., 
Kravtsov et al. 2004) or by an AGN (Croton et al. 2006). 
However, even large variations of these parameters have 
not been successful in solving the problems significantly 
for all galaxy types. Another kind of solution is resorting 
to the modified Newtonian dynamics (MOND), proposed 
by Milgrom (1983) as an empirical modification of grav- 
ity, when the generated acceleration falls below a universal 
value Qo ~ 2 x 10~ 10 m s~ 2 . In this model, there is no DM 
anymore, but the visible mass in the inner parts of galaxies 
produces a much boosted gravity force in the outer parts, 
with a longer range effect. Bekenstein & Milgrom (1984) 
developed a self-consistent Lagrangian theory, where the 
Poisson equation is transformed into: 



V[/4|V$|/a )V$] = AwGp, 



(1) 



MOND regime. Far in this regime, and assuming some sym- 
metry (spherical, cylindrical, or plane) it can be shown that 
the MOND acceleration qm satisfies the relation (Brada & 
Milgrom 1995): 



9 m = a o9N, 



(2) 



where /i(x) is a function that is equal to unity at large x 
(Newtonian regime), and tends to x when x << 1 in the 



where gjy is the Newtonian acceleration. This model has 
large success at galactic scale, in particular explaining all 
rotation curves of galaxies, and naturally the Tully-Fisher 
relation, as developed in the excellent review by Sanders & 
McGaugh (2002). 

Interest has grown in the MOND theory since the propo- 
sition by Bekenstein (2004) of a Lorentz-covariant theory 
(TeVeS), able to replace general relativity, accounting for 
gravitational lensing and passing elementary tests of grav- 
ity in the solar system. Simulations have been attempted to 
explore the large-scale structure formation, with encourag- 
ing results (Knebe & Gibson 2004; Nusser & Pointecouteau 
2006). More recently, weak lensing observations of the bul- 
let merging cluster 1E0657-56 (Clowe et al. 2006) claim that 
the spatial separation between the main baryonic compo- 
nent (X-ray gas) and the total mass shows direct evidence 
for the existence of collisionless DM. They find that any 
modified gravity model, considering only the baryonic mass, 
fails to reproduce the observations. However, Angus et al. 
(2006) have re-analyzed these observations in the context of 
modified gravity and show that the data are also compat- 
ible with the Bekenstein model of MOND, in which some 
collisionless dark matter exists under the form of ordinary 
hot neutrinos of 2 eV. 

The most stringent constraints on the choice of the in- 
terpolation function (x(x) are expected to be obtained on a 
small scale however. To better fit the rotation curve of the 
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Milky Way, the function p(x) = x/(l + x) has been pro- 
posed by Famaey & Binney (2005) , in place of the empirical 
initial function p(x) = x/(l + a; 2 ) 1 / 2 . In addition, physical 
constraints and the external field effect further reduce the 
choice of the interpolating function (Zhao & Famaey 2006). 

Since the motivation of MOND and its best success 
concern the galactic scales, and in particular the rotation 
curves fit without dark matter, more physical constraints 
should be explored at these scales. In particular, the sta- 
bility of spiral galaxies in this model, the secular evolution 
taking into account spiral waves and bars have to be in- 
vestigated, to compare the dynamical behavior of a typi- 
cal galaxy in the Newtonian CDM model and the MOND 
frame. Brada & Milgrom (1999, hereafter BM99) have be- 
gun to tackle this problem, and have shown that the Toomrc 
Q-parameter could be chosen lower than in the Newtonian 
case, to obtain the same stability level. The modified accel- 
eration provided a comparable stability level with respect 
to bars as does a dark matter halo in the Newtonian case. 
There are, however, limitations in their model, since they 
considered infinitely thin discs and ignored the z-structure, 
acceleration, and dynamics, which are very different in 
Newtonian and MOND regimes. 

In this work, we present numerical simulations of sev- 
eral spiral galaxy models, representing the whole Hubble 
sequence and a large mass range, in both CDM Newtonian 
and MOND models. The goal is to find specific tests and 
constraints to the gravity theory, to be applied on a global 
statistical basis and confront them to the observations. The 
diagnostics are to be found in the bar frequency, the spiral 
morphology, the thickness of discs and their box/peanut 
shapes, the surface density profiles, and the angular mo- 
mentum distribution. In this first approach, pure stellar 
discs are considered, while gas and star formation will be 
investigated in a future work. In the next section, we de- 
scribed the numerical code developed to solve the difficult 
problem of MOND dynamics, and in Sect. 3 the analy- 
sis and diagnostics we applied to the simulations results. 
Initial conditions for spiral galaxies described in Sect. 4, are 
selected to be as close as possible in the plane for the two 
compared models: in particular they have the same radial 
baryonic distribution and the same rotation curve and ve- 
locity dispersion. Results are presented in Sect. 5 and then 
discussed in Sect. 6 to emphasize the fundamental differ- 
ences in galaxy evolution for the two competing dynamics. 



2. Numerical model 

The non-linearity of the MOND gravity leads us to use 
different techniques than the usual ones for the potential 
solver (or force solver). 

2.1. Multigrid (MG) potential solver 

The modified Poisson equation is a non-linear elliptic par- 
tial differential equation (PDE). This kind of equation can 
be solved efficiently using multigrid (MG) techniques. We 
have written an N-body code in which we implemented 
a full multigrid algorithm (FMG) with full approxima- 
tion scheme (FAS) for the potential solver (see Numerical 
Recipes, Press et al., 1992). Brada and Milgrom (BM99) 
used such a code to solve ([!]). 
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Fig. 1. Full multigrid (FMG) algorithm is used to accel- 
erate the convergence in the resolution of the modified 
Poisson equation (see text). 



Up to some point, the code works like a particle- 
mesh(PM) code. Particles evolve in a 3D Cartesian grid. 
Density is computed using the cloud in cell interpolation, 
the potential is deduced by MG techniques, the equation of 
motion is solved by the leapfrog scheme. The only difference 
from a classical PM-code occurs in the potential solver. 

The MG computes the solution on finer and finer grids 
(Fig. [T]) by calculating correction terms on each level and 
converges even more quickly than by solving the same equa- 
tion directly on the finest grid. We use the Gauss-Seidel 
relaxation with red and black ordering (Press et al. 1992) 
to solve the system of equations obtained by discretisation. 
This step is called smoothing. To go from the grid level n 
to n + 1 we make a tri-linear interpolation (prolongation 
operator, P), and inversely, the full- weighting operator (R) 
is used to go to the level n + 1 to n. The number of pre- 



/post-relaxations were chosen to v pre 
Here is the discrete form of ([1]): 



2 and u 



post 



1. 



47rG/9 ijifc = (3) 

{4>i+l,j,k — </>i,j,fc)MMi — (4>i,j,k — </>i-l,j,fc)M£i 
+ W,j+l,k ~ 4>i,j,k)PM 2 — {4>i,j,k — 4>i,j-l,k)PL 2 



+ (4>i,j,k+l ~ <f>ij,k)PM 3 



with Pij t k and <pij t k the spatial density and potential dis- 
cretized on a grid of step h, PM n and /if,,, the value of 
p{x) at points M; and Li (Fig. The gradient com- 
ponent (d/dx,d/dy,d/dz), in fx(x), are approximated by 

, <t>(B)-<KA) 4>(i)+HH)-4>(K)-4>(J) MC)+<t>(D)-ME)-<t>(F) \ 

y h ' ih ' ih )i 

it is the stable numerical scheme proposed in BM99. 
In the Newtonian case, the interpolation function is con- 
stant: p(x) — 1, so that Eq. becomes: 

4wGpij t k = (4>i+i ,j,k + 4>i-lJ,k + <t>ij+l,k + <t>i,j-l t k (4) 
+ <f>i,j,k+l + <k,j,k-l - 60ij,fc)//l 2 

We recognize the discrete Poisson equation with a 7-point 
Laplacian stencil. 

To avoid the usual 2-body relaxation in simulated galax- 
ies with insufficent number of particlees, the gravitational 
potential is softened through a convolution with a gaus- 
sian function (er = 1.2 cells). This value of the softening 
suppresses efficiently high spatial frequency noise, without 
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Fig. 2. Discretisation scheme of the modified Poisson equa- 
tion proposed in BM99. Density and potential is calculated 
on the grid nodes. The gradient components in n(x) are 
estimated at the Li and Mj points. 



introducing any bias (Dehnen 2001; Zhan 2006). For the 
following simulations, the calculation was made on a 257 3 
grid. The radius of the simulation box is 50 kpc, generally 
the galactic disc is truncated at 20 kpc. The same code can 
solve Poisson and Modified Poisson equations (it is just the 
coefficients of the PDE that are constant in Newton and 
variable in MOND). 



2,2. Boundary conditions 

In classical PM-code, fast Fourier transform (FFT) implies 
periodic boundary conditions. For isolated galaxy simula- 
tions, the interactions with periodic images can be sup- 
pressed using screening masses (James 1977). What kind 
of boundary conditions are possible to use for simulations 
of isolated galaxies with MG codes? This point is not trivial 
and is not developed in BM99 code. It is particularly im- 
portant in MOND where the gravitational potential scales 
as log(r) far from the galaxy. It might appear similar in 
Newton gravity with the dark matter halo, but the latter 
is nearly spherical in general, and the influence of mass 
exterior to the box is considered negligible. In MOND in 
the contrary, the potential at large distance is due to the 
baryonic disc (with spiral arms or bar). 

Periodic conditions are not realistic at this scale and 
using a box eight time larger to suppress the periodic im- 
ages is too costly in CPU time. The most natural way is 
to use isolated boundary conditions. But this supposes we 
know the potential at the border of the box. To solve this 
problem we have to make an approximation. We use the 
MOND formula (Eq.© in the deep MOND regime), which 
links the MOND acceleration to the Newtonian accelera- 
tion. If fi(x) = x/(l + a;), one has more generally: 



°MOND - a Newt^ ( a Newt/ a o) 
with 

v{x) = 0.5 (l + v 7 ! + 4/x) . 



(5) 



(6) 
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Fig. 3. Plot of | V(j>newt | versus (f> n ewt for a barred 
galaxy with grand design spiral. The crosses represent 
I ^4>newt \ = <Pnewt/GM (this is the spherical approximation 
where the galaxy is approximated by a point mass). 



It is critical to use this expression directly in the Newtonian 
code principally because it is not true if the system has no 
symmetry (planar, cylindrical, or spherical). It is a concern 
even for an isolated galaxy. During its evolution, spiral arms 
and bars are formed and destroyed. The particle configura- 
tion is then not symmetric. However, we are interested only 
in the outer parts of the galaxy. We just need to determine 
the MOND potential on the boundary of the simulation 
box, that is far from the galactic center and its gravita- 
tional instabilities. 

Brada and Milgrom (1995) proposed a test to check if 
the approximation ajyjQND — /( a newt) i s justified. They 
showed that | V(f> n ewt \ must be a function of (j> n ewt out 
of the disc. So by plotting | V<j) new t | versus 4> n ewt at dif- 
ferent positions in the box simulation we obtain Fig. [3l 
The approximation is good for low <f> that is far from the 
galactic center. This is expected since the potential is more 
spherical. Hence, our solution to solve the boundary con- 
ditions problem is to compute the Newtonian potential by 
the FFT technique on a larger grid, then the Newtonian 
acceleration on each edge of the simulation box. Finally we 
use the MOND formula to deduce the MOND acceleration 
and compute the MOND potential on the border. 

In this way, we obtain boundary conditions that are not 
fixed in time and that are not required to be homogeneous. 
We can take a small perturbation to the spherical symmetry 
like the disc or bar shape of the galaxy into account. Even if 
the correction is not very important to dynamical evolution, 
this makes the code more realistic. 



2.3. Tests 

We have made several tests to check the validity of the solu- 
tion obtained by the MG technique: the analytical solution 
of a mass point, the Kuzmin disc. We present here a more 
demanding test for a totally non-symmetric system. It is 
the potential of a galaxy where a bar is formed during the 
simulation. We compute the MOND potential on the one 
hand with the MG technique and on the other hand with 
a classical relaxation scheme (NAG library). The second 
method is very inefficient (several hours when it takes less 
one minute with MG). We plot the potential along the bar 
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Fig. 4. Test of the MG solver compared to the NAG routine 
DOECF to solve modified Poisson equation for a barred 
system with a grand design spiral structure; (1) potential 
perpendicular to the bar, (2) potential parallel to the bar. 

and perpendicular to the bar. The two different methods 
are in complete agreement. This plot also demonstrates the 
high symmetry at the outer boundary of the galaxy. Even 
if there is a bar (bkpc), the potential at 25 kpc is quasi- 
spherical (within a few %). 

We have tested the dynamical evolution of a stellar disc 
in Newtonian gravity with the MG code compared to a 
classical FFT-PM-code. We obtained the same result for 
the time evolution of the bar. The code has been paral- 
lelized in open-MP since the G-S relaxation with red and 
black ordering allowed this. All red cells can be updated 
independently and this the same for the black cells. 

3. Analysis techniques 

3.1. Fourier analysis in the galactic plane 

The potential in the galactic plane is developed with the ba- 
sis of the cosine function (<j> m (r)) and a phase term ((/> m (r)): 

$(r, 6) = $ (r) + $ ™( r ) cos i md ~ <f>m{r)} 

m— 1 

to calculate the maximum strength (Q m ) of the m mode. 
We use the maximum force ratio: 



Q, 



max 



with the radial force: 



Fr 



dr 



and the tangential force: 



max(Fg. 



1<9$ 

= 1 77w 

— $m(r). 
r 



The bar strength is the maximum strength of the mode 
m = 2. In general, the phase term, fa, gives the rotation 
speed of the bar flf, and the derivative: 



But the mode m = 2 could correspond to a two-arm spi- 
ral structure. Then, it is more informative to calculate the 
Fourier transform fair, O) from fa(r,t). One can distin- 
guish the angular velocity (f2) of a structure like a bar or 
a spiral arm versus the radius (r). A bar is identified by a 
solid rotation in the central part of the galaxy: 

fi(r) = est = fib 



3.2. Resonance 

We estimate the position of resonant orbits using the 
epicyclic approximation (Fig. [TO]) . To do that, we need to 
determine the angular velocity of the stellar disc (11), 



the epicyclic frequency (k), 



1 rf$ 
r dr 



AQ 2 



2 dn 2 

k = r— — 
dr 

and the vertical frequency (f z ), 
2 _ d 2 $ 



3.3. Heating 

The heating of disc is computed by averaging the radial 
velocity dispersion, cr r (t), normalized by the initial a cr u 
(see Sect. 14. 2p inside the 5 kpc of the galaxy, giving the 
averaging Toomre parameter: 

Qt = ( )r<5 kpc 

O crit 

with a cr it, the critical velocity dispersion derived from the 
Toomre stability criterion, 



&crit 



3.36GS 



dt 



E is the stellar surface density. 



3.4. Units 

In our code we use a unit system where the universal 
constant of gravity is: G = 1, and the mass unity is 
U m = 2.26 10 9 M Q . The length unity is U r = 1.02 kpc 
and the velocity unit is U v = 100 fcm.s -1 . The time unit 
is Ut = 10 Myr. In this paper, when the units are not 
indicated, they are in this unit system. 



4. Initial conditions 

To study the stability differences for galaxies in the MOND 
and DM models, we construct a sequence of galaxies from 
early type to late type (Sa, Sb, Sc, Sd). Each type of galaxy 
corresponds to a set of two model galaxies, one for MOND 
and the other for DM. A galaxy, for a given type, has the 
same spatial density for the baryonic disc and bulge in the 
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Run 


M d 


M b 


ad 


M h 


b h 


Sa 


40 


12.65 


4 


206.4 


14.8 


Sb 


30 


5 


5 


173.7 


14.6 


Sc 


20 


2 


6 


148.8 


14.5 


Sd 


to 




6 


129.7 


13.1 



Table 1. Parameters. The characteristic length of the bulge 
for Sa, Sb, Sc galaxies is 1 kpc. For the characteristic height 
of the Miyamoto-Nagai disc, we choose bd/ad = 1/10. The 
given mass is the truncated mass. The Toomre parameter 
value is the same in the DM and MOND model: Qt = 2. 
The disc is made of 2.10 5 particles, the mass of the bulge 
particles is equal to the mass of the disc particle. The mass 
of dark matter particles is three times the mass of the disc 
particle. 



1 

0.1 
0.01 
0.001 
le-04 
le-05 
le-06 



Newton 
MOND 
Newton+DM 



two models. The stellar disc is modeled by a Miyamoto- 
Nagai disc: 



b 2 M d \ a d R 2 + (a d + 3^T¥ d )(a d + + 



-in 



R 2 + {a d + ^z 2 + b 2 ) 2 (z 2 + 62)3/ 



5/2 



with Md the mass of the disc (at infinity), ad and bd 
the characteristic length and height, and the bulge by a 
Plummcr sphere: 



Pb 



-5/2 



with Mb the mass of the bulge (at infinity) and b^ the char- 
acteristic length. 

4.1. Rotation curves 

From an observer's point of view, a galaxy must have the 
same rotation curve in MOND and in DM. The shape of the 
rotation curve is imposed by the MOND model. To obtain 
the same rotation curve (in the galactic plane) in DM we 
adjust a Plummer dark matter (live) spherical halo to fit 
the MOND rotation curve. Parameters of the dark matter 
halos are given in Table [TJ The error on the fit parameters 
is about 2 %. 

4.2. Velocity dispersion 

We used the same value for the Toomre parameter (Qt) for 
MOND and DM. The radial velocity dispersion is initialized 
by: 

The tangential velocity dispersion (erg) is deduced from the 
epicyclic approximation, 



(70 = C, 



2VL' 



For the vertical velocity dispersion (cr z ), the hydro- 
static equilibrium of an isothermal infinite stellar disc in 
Newtonian gravity gives: 

a 2 z = HirG^r) 

where H is the characteristic height and £(r) is the surface 
density. 



Fig. 5. Vertical structure at the gravitational equilibrium 
of a disc in Newtonian gravity, MOND, and Newtonian with 
a dark matter halo. 



We have calculated the vertical density profile (p(z)) in 
MOND for an isothermal infinite stellar disc (we consider 
the problem in one dimension). The equivalent pressure of 
the gas of star is P — pa 2 . The gravitational potential <& is 
given by the modified Poisson equation. The gravitational 
equilibrium (VP = — /9V$) is obtained when: 



d 

dz 

with 
d$ _ 

dz 



d$/dz \ \ d$ 
dz 



= 4nGp 



1 dp 

p dz 



(7) 



(8) 



a z is a constant of z but varies with r. We solved numer- 
ically Eq.0 the result is plotted in Fig.O It shows also the 
vertical profile in Newtonian gravity: p(z) — po sech 2 (z/H), 
and in Newtonian gravity with a dark matter halo. For 
this plot, we selected po = 2.3 10 -6 Mqy.kpc" 1 , which is a 
typical value of the outer disc. The dark matter halo is a 
Plummer sphere with Mh = 6 IO^Mq and bh = 15 kpc 
(r = 8 kpc). In our model we choose H = est in the DM 
model as well as in the MOND model. We want to keep the 
same initial height for a galaxy in DM an MOND. Figure 
[5] shows that the vertical profile in MOND or in Newton 
with a dark matter halo are quite similar. The initial verti- 
cal velocity dispersion in MOND is the same as in the DM 
model. The stellar rotational velocity is not exactly the cir- 
cular velocity (v c ), but v c — v a where v a is the asymmetric 
drift deduced from Jeans equations applied to an infinitely 
thin disc. The system is relaxed initially in its axisymmetric 
potential to have a well stable virialized initial state. 

5. Results 

5.1. Bar growth 

5.1.1. Dark matter model 

Figure [6] (left panel) shows the evolution of an Sa galaxy 
in the DM model. At first sight, the initial Miyamoto- 
Nagai disc develops a bar instability in several Gyr. The 
bar length grows to 6 kpc until 2 Gyr (Fig. [8]); its shape is 
rather squarish. For this run, we do not clearly see a grand 
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Fig. 6. Bar growth of Sa type (Qt = 2) in the DM model (left panel) and the MOND model (right panel). The size of 
the box is 80 kpc x 80 kpc. In the DM model, the bar develops in several Gyr. It can be noticed that the bar is surrounded 
by a ring at the end of the simulation. Particles are confined in the disc. In the MOND model, the bar appears quite 
rapidly (in less than 1 Gyr) , and a lot of particles are spread out around it up to 30 kpc. 



0.25 



0.15 





0.05 



3 4 5 
Time (Gyr) 

Fig. 7. Time evolution of Q m for m = 1, 2, 3, 4, 8 of galaxy Sa in DM model (left) and in MOND model (right). In the 
DM model, the bar strength increases progressively compared to the MOND model where the bar reaches its maximum 
after 1 Gyr. The same drop appears at t = 2.5 Gyr in DM and t = 4.5 Gyr in MOND. After that, the bar strength 
increases again in the DM model, but not in MOND. 



DM bar length 
MOND bar length 

DM corotation 
MOND corotation 




3 4 5 
Time (Gyr) 



Fig. 8. Bar length and corotation radius in the DM and 
MOND models. The bar length is defined by the radius 
where the bar strength (Q2M) is equal to half its max- 
imum. The difference between the end of the bar and 
the corotation radius in MOND is constant and smaller 
than in the DM model. The corotation radius in the DM 
model is shifted outward because the bar slows down 
(see Sect.lQ]). 
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design spiral structure during the bar growth. They exist, 
but they are more visible for a colder disc (Qt = 1-5). 
After t = 2 Gyr, transient spirals are developed between 
10 kpc and 20 kpc, while a ring appears at the end of the 
bar. The bar length continues to grow until 4 Gyr, as the 
ring extends, too (6 — 7 kpc). During the period between 
5 Gyr and 8 Gyr, the bar changes its morphology and takes 
a butterfly shape. Its length does not increase contrary to 
the ring that has extended to 10 kpc. Spirals that devel- 
oped at t = 2 Gyr have driven particles to a pseudo ring at 
15 kpc, for t = 4 Gyr, to 25 kpc, at t = 8 Gyr. The nature 
of these rings will be discussed in Sect. I5.2T21 

Figure[7](left) displays the maximum strength of the bar 
as a function of time. One can distinguish three parts on 
this plot. First, Q2 begins to increase until 2 Gyr. Then, this 
growth stops suddenly, the bar strength drops in 500 Myr. 
After 2.5 Gyr the bar strength grows again until 5 Gyr and 
appears to be constant until 8 Gyr. 

5.1.2. MOND model 

In MOND, the same galaxy (with the same value for Qt) 
shows quite different structures (Fig. right). First, a 
multi-spiral pattern appears rapidly after 0.4 Gyr to give 
place to a grand design two-arm spiral at 0.6 Gyr. This 
spiral structure persists until 3 Gyr. The spiral arms have 
spread out particles up to 30 kpc. After 4 Gyr, the bar 
begins to be rounder and weakened. No rings are clearly 
visible in the MOND simulations. 

Qi(t) is plotted on Fig. [7| (right), the bar develops very 
soon (1 Gyr), compared to the DM model. The bar strength 
is constant until 4.5 Gyr where a drop occurs suddenly (like 
in DM at 2.5 Gyr). Afterwards, the bar strength remains 
low until the end of the simulation. However, the bar length 
is still constant (6 kpc) even if the bar strength is weakened 
(Fig.®. 

5.2. Pattern speed and resonance 

5.2.1. Pattern speed 

The bar pattern speed is represented on Fig. [9j still for 
the Sa-type galaxy (Qt = 2), for the MOND model, and 
for the DM model with a live and analytic halo. The bar in 
DM with a live halo is considerably slowed down during the 
simulation (25 km.s~ l .kpc -1 to 10 km.s~ l .kpc~ l ), while in 
MOND, the pattern speed is constant (25 fcm.s -1 .kpc -1 ). 
This plot emphasizes the dynamical friction effects experi- 
enced by the stellar bar against the DM particles. To con- 
firm this result we perform a second simulation with the 
DM model using an analytical dark matter halo instead of 
a live halo. In this case, the pattern speed of the bar is still 
constant and corresponds to the MOND result. 

5.2.2. Corotation & ILR 

The bar pattern speed determines the position of resonant 
orbits in the reference frame of the bar rotating at f2f,. 
Because of velocity dispersion, stars do not just have a cir- 
cular motion around the galactic center, they oscillate with 
the epicyclic frequency k (parallel to the galactic plane); 
likewise, they oscillate in z with the frequency v z . In most 
numerical simulations (with gas and dark matter) and in 
galaxies where it was possible to determine the bar pat- 
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Fig. 9. Rotation speed of the bar for MOND and DM mod- 
els. In MOND, the bar turns with the same velocity from 
the beginning to the end. In DM, particles of the halo slow 
down the rotation of the bar by dynamical friction. 

tern speed, the bar extends to its corotation (e.g., Buta & 
Combes 1996). 

In the DM model, during the two first Gyr, the bar pat- 
tern speed is about 25 km.s~ l .kpc~ l so that the corotation 
is nearly 13 kpc (Fig. OH top), while the bar ends at 7 kpc, 
but as the bar slows down, the corotation is shifted out to 
20 - 25 kpc at t = 8 Gyr (Fig. [TQl middle). Hence the ring 
surrounding the bar from t = 2 Gyr until the end of the 
simulation does not correspond to the corotation resonance 
as it might be expected. The epicyclic approximation (Fig. 
[TU1 middle) indicates that we should have an outer and an 
inner Lindblad resonance (OILR, IILR), that is where Of, 
(the bar pattern speed) intercepts f2 ± k/2. The IILR is 
located very near the center of the galaxy. The OILR is 
about 12 kpc, where the ring is observed at the end of the 
simulation. The bar appears to end nearly at the OILR. 

Between the two ILR X2, orbits must exist and destroy 
the bar. Figure [IT] displays the velocity field in the refer- 
ence frame rotating with the bar, and the potential outline 
indicating the bar orientation. Corotation is well identified 
when vectors change orientation (r ~ 25 kpc). The trajec- 
tory of particles are parallel to the bar potential like x\ or- 
bits, not perpendicular (like X2 orbits) . We have performed 
an orbit analysis to determine the existence of x\ and X2 
orbits in this bar potential. The result is that the bar is not 
dominated yet by X2 orbits. We have launched particles in 
the bar potential rotating at fit, = 10 fcm.s -1 .kpc -1 . The 
value of the Jacobi's integral, 

Ej = h 2 + $ - ~ I fl h x r | 2 

of particles varies between h rnin = —38 (in our system 
units, G = 1 see Sect. I5.2.2|) . the bottom of the poten- 
tial well, to h ma x = — 13, the potential nearly the corota- 
tion. X2 orbits exist between r = 3 — 3.5 kpc in an en- 
ergy range about —20 < Ej < —18. For lower like 
5 fcm.s -1 . kpc -1 , X2 orbits appears clearly and more fre- 
quently between r = 2.5 — 9 kpc in an energy range of 
—26 < Ej < —13. Figure [TO] gives just an indication on the 
resonance with the epicyclic approximation. It can be noted 
that the drop in bar strength between 2 Gyr and 2.5 Gyr 
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is correlated with ILR formation and its analogue in the 
z-direction (peanut formation, see the next section). In the 
MOND model the bar always ends near its corotation (Fig. 
[8]), while the DM bar is relatively shorter. 

5.2.3. Vertical resonance and warp 

It can be shown that v z and k have a similar evolution. An 
equivalent resonance of the ILR exists in the z-direction if 
fib = fl — v z j2. When particles resonate both in the plane 
and perpendicular to it, their vertical oscillations can be 
amplified and a peanut shape results. 

The drop in the evolution of the bar strength coincides 
with the peanut resonance, in DM as well as in MOND. 
At this moment particles are elevated out of the galactic 
plane. Stars are less bound and orbits become more oval, 
the bar strength is thus weakened. Figure [T2l illustrates the 
moment when the peanut is formed. At t = 2.5 Gyr in 
the DM model, and t = 4.5 in the MOND model, particles 
between 2 kpc and 8 kpc resonate and are detached from 
the galactic plane. 

To confirm that the drop in the bar strength is really 
due to the peanuts, we have performed a 2D-simulation for 
the DM and MOND models. In this case the bar strength 
is not weakened during its evolution Fig. 1131 There is no 
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Fig. 11. Velocity field in the reference frame rotating with 
the bar. Vectors are parallel to the isopotential lines, indi- 
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Fig. 13. 2D simulation of the Sa galaxy in DM and MOND 
models. The bar strength increases until a maximum and 
remains roughly constant. 



drop at t = 2.5 Gyr in the DM model or at t = 4.5 Gyr in 
the MOND model. 

Like ILR and corotation in the DM model, the position 
of the z-resonance is also shifted out when the bar pat- 
tern speed decreases. In our simulation, dynamical friction 
acts quite progressively on the bar. The radial shifting of 
the peanut lobes is continuous. Martinez- Valpuesta et al. 
(2006) have observed similar phenomena in the formation 
of a peanut galaxy, but they distinguish two episodes. They 
obtained a short bar due to chaotic orbits that appear be- 
tween the ILR and vertical-ILR. 

Figure [TJ] shows several edge-on views of the galaxy 
in the DM and MOND models. In MOND simulation, 
the galactic disc is more easily warped than in DM. The 
disc begins to take a U-shape, to finally flare. The ratio 
h/h r , where h is the equivalent characteristic height and h 
the characteristic length for an exponential disc, is about 
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Fig. 12. Peanut formation in the DM (top) and MOND (bottom) models. These plots display {z 2 ) 1 / 2 of particles at 
several radii in a function of time. The peanut appears at t = 2.5 in the DM model and t — 4.5 in the MOND model. 
At these times, particles get out of the galactic plane in the range between 2 — 8 kpc. As the bar slows down in the dark 
matter halo, more and more particles resonate, and the peanut lobes extend to 12 kpc at t = 8 Gyr. In the MOND model 
the peanut stays at the same place all along the simulation. 





Fig. 14. Edge-on view of an Sa galaxy, shows the charac- 
teristic peanut shape. In the DM model (left), the position 
of the peanut lobes are radially shifted out as the bar slows 
down. In the MOND model (right), the peanut keeps the 
same size, one can notice the warp and flaring of the disc. 



0.26 in the MOND model and 0.22 in the DM model at 
r = 25 kpc. It is not as different as can be expected be- 
cause h r in MOND is larger (7 kpc) than in the DM model 
(5.5 kpc). Particles are ejected radially further than in the 
DM model because of the angular momentum transfer (see 
Sect. 15 .4[) ; the disc is thus less compact. The origin of the 
flare comes from the vertical velocity dispersion that is 
more important in MOND for outer regions than in the 
DM model (Fig. [T5]). In MOND, vertical instabilities are 
developed because of the self-gravity, and the disc heats 
more. 

We have seen that the weakening of the bar coincided 
with the peanut's occurrence. Let us note that the peanut is 
not the only way to weaken a bar in a pure stellar disc. If the 
disc is too cold, it develops a bar instability very quickly 




Fig. 15. Vertical velocity dispersion (<r z ) of the disc for the 
DM and MOND models at t = 8 Gyr. In DM, the peanuts 
makes a z increase around 12 kpc. After 15 kpc, a z in MOND 
is larger than in DM. 



so that the stars have no time to follow a typical orbit 
with a constant bar strength. The corresponding disordered 
motions of the stars weaken the bar. This can be shown in 
a 2D (planar geometry) simulation (Fig. [T6|) . 

5.3. Heating 

The problem now is to understand why the bar strength 
continues to increase after the z-resonance in DM and not 
in MOND. Part of the explanation can be found by fol- 
lowing the evolution of Qt for the two models. The value 
of the Toomre coefficient indicates the heating rate of the 
disc. In these simulations, Qt starts at a value of 2 in the 
whole disc. The evolution in the DM model and MOND 
model (Fig. [T7|) is differentiated from the beginning like 
the evolution of the bar strength. 



Heating in the DM model. In the DM model, the disc heats 
progressively. When the peanut is forming, it weakens the 
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Fig. 16. 2D simulation of a cold stellar disc (Qt = 1) evolv- 
ing in a small dark matter halo. The bar is formed quickly 
in a few galactic rotation. Bar weakens because of non- 
adiabatic bar growth (see text). 

bar and Qt ~ 2.7. This value is not enough to avoid bar 
formation, which is why bar strength increases again. A 
disc in DM model needs a value of about 4 for Qt to be 
stable and not form a bar. 

Heating in the MOND model. In MOND, Q T increases to 
3.5 in a few Gyr. The apparition of a z resonance weakens 
the bar strength. At this time, the MOND disc is thus more 
stable because particles have more velocity dispersion. The 
bar strength does not increase anymore. We have performed 
another MOND simulation with Qt = 3.5 from the begin- 
ning. A weak bar is formed with a strength of about 0.12. 
That corresponds with the bar obtained at the end of the 
simulation with Qt = 2 initially. 

In MOND, all the matter participates to the dynam- 
ics; the galactic disc is completely self-gravitating, hence it 
heats up more than a disc rotating in a dark matter halo. 
The fact that no ring is clearly visible in the MOND simu- 
lation can be understood since the disc is hotter and might 
not sustain these features. The pitch angle of a density wave 
depends on Qt- For a hot system, the theory predicts that 
the spiral will be more open than for a cold one; it is thus 
more difficult to form a ring. 

5.4. Angular momentum exchange 

Another crucial point for the bar formation is the exchange 
of angular momentum. For the bar to grow, particles of the 
disc have to lose angular momentum to fall in the inner 
region and have an elliptical orbit instead of a circular one. 
Angular momentum can be exchanged between the inner 
and outer parts. 

Angular momentum and dark matter. In DM, the halo can 
receive angular momentum from the disc. It is well illus- 
trated in Fig. [THJ the disc loses about 30% of angular mo- 
mentum in the halo. In other terms if the halo increases its 
angular momentum, it will be less compact and will inflate. 
Figure[Tn]represents the time evolution of several radii com- 
prising a fraction of the mass between 10% to 90% by 10%. 
One can notice the expansion of the radius below 60%. At 
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Fig. 17. Evolution of the Toomre coefficient value in func- 
tion of time. The disc heats up in a few Gyr in the MOND 
model because of a complete self-gravity. 
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Fig. 18. Vertical component of the angular momentum in 
DM model. L z is exchanged from the disc to the halo during 
the bar formation. A non negligible proportion of particle 
(dark matter particles) escaped from the simulation box. 
Their motions are treated in a Keplerien potential as if all 
the matter inside the box is concentrated into a point mass. 



t = 1 Gyr, 90% of the mass is included in a sphere of 29 kpc; 
at t = 8 Gyr this mass is in a sphere of 32 kpc radius. 

This exchange, which is efficient after 2 Gyr, contributes 
to bar formation (especially after the drop). Then the disc 
is not too hot and density waves can propagate. One can 
notice that the core of the halo is unaffected. 

Angular momentum in MOND. In MOND, the disc does not 
lose lot of angular momentum. Angular momentum is ex- 
changed inside the disc itself. Figure [2D] shows the evolution 
of the same radius (10 — 90%) of the mass versus time for 
the disc. The inner part of the disc loses angular momen- 
tum as expected because of bar formation (contraction of 
the disc below 8 kpc), and the outer part of the disc re- 
ceives angular momentum from the inner region (90% of 
the mass is inside a 15 kpc at the beginning and extends at 
20 kpc at the end). This occurs during the first 3 Gyr, and 
the transfer is mediated by the spiral arms seen in Fig. [SJ 
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Fig. 19. In the DM model, the dark matter halo inflates 
because of angular momentum exchange from the disc to 
the halo. 
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Fig. 22. Bar pattern speed in DM (left) and MOND (right) 
models. Dynamical friction effects in the DM model slow 
down the bar. In MOND pattern speed is still constant. 
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Fig. 20. In MOND, the outer part of the disc inflates while 
the center contracts. Angular momentum is transferred be- 
tween these two regions during the bar growth. 



They evacuate angular momentum from the inner part to 
the outer part of the disc, and spread out particles around 
the disc. This is possible when the disc is not too hot. After 
this phase there is a saturation when the disc becomes sta- 
ble and no density wave can propagate and increase the bar 
strength. 

5.5. Dark matter compared to MOND along the Hubble 
sequence 

A series of simulations have been run to explore the pa- 
rameters of galaxies along the Hubble sequence according 
to Table. [T] In the DM model, from the early-type to late- 
type galaxies the ratio between the visible mass and the 
dark matter inside the optical radius increases. Thus, late- 
type galaxies are less self-gravitating than early types, so 
they are more stable (Fig. [2Tj) . 

In MOND, galactic discs are cold and form a bar in a 
few Gyr whatever their type. The evolution scheme seen 
for the Sa type is reproduced for the Sb, Sc, and Sd type 
too. Even if the disc in MOND is cold and unstable at 



the beginning, it heats quickly and stabilizes itself along its 
evolution fFig. l2Tj) . 

Peanuts are formed for galaxies with a sufficiently mas- 
sive bulge like Sa and Sb galaxies. In this case, the peanut 
occurrence weakens the bar. But in the DM model, the bar 
strength increases again because of the angular momentum 
transfer between the disc and the dark matter halo. While 
in MOND, the bar strength keeps low, and the disc heats 
up because of instability and stabilizes itself. 

For Sc and Sd galaxies in the MOND model, the bar is 
formed too quickly (a few galactic rotations) because the 
system is too unstable. The stars have no time to settle in 
orbits supporting the bar at a given bar strength, since the 
orbital structure of the bar varies on a time scale shorter 
than the orbital period. These galaxies present a strong bar 
during a short time at the beginning of simulation to finish 
with a weak bar. 

The pattern speed of the bar is plotted in Fig. [22j In 
MOND the pattern speed is always constant for a given 
galaxy. Early-type galaxies have a higher bar pattern speed 
than late types (the disc is more massive). In DM, the bar is 
always slowed down by dynamical friction due to the halo. 
Late-type galaxies need more time to form a bar (Fig. |2"T|) , 
so their pattern speeds are less slowed down than for early 
types. 

We have made a statistical study of the bar strength 
for typical galaxies of the Hubble sequence. Figure [23] shows 
the bar frequency obtained using the time spent by a galaxy 
with a given bar strength. Two tendencies are clear from 
this plot. First, galaxies in the MOND model have stronger 
bars (Q 2 > 0.25) than in the DM model. The MOND 
discs are more unstable at the beginning so they form a 
strong bar very quickly. Secondly, there is a hole at low bar 
strength in the MOND model that is not present in the DM 
model. This is due to the dark matter halo that stabilizes 
the disc at the beginning, so it takes more time to a galaxy 
to form a bar. The bar strength distribution obtained from 
the observations presents some characteristics that are re- 
produced with the MOND model. In particular, there is a 
small proportion of galaxies with a very weak bar, and a 
few galaxies have very strong bar (e.g., Block et al. 2002 
and Whyte et al. 2002). 
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Fig. 21. Evolution of the bar strength in the DM (left) and MOND (right) models. Late-type galaxies in the DM model 
are more stable than early types; they need a larger proportion of dark matter to obtain the same rotation curve as in 
MOND. Late-type galaxies in MOND are too unstable; the bar destroys itself. 
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Fig. 23. Bar frequency in the simulated Hubble sequence 
in MOND (top) and in DM (bottom). Bars are stronger 
with MOND and there is a dearth of galaxies without bars 
in MOND, but not in the DM model. 



6. Discussion and conclusion 

In this paper, the dynamical evolution of pure stellar discs 
in MOND is compared to Newtonian gravity with DM, us- 
ing numerical simulations. We have developed an N-body 
code that solves the modified Poisson equation in three di- 
mensions using MG technique for the potential solver. The 
simulations in the DM models have been performed with 
the same code by solving Poisson equation with the same 
MG technique. 

For isolated galaxy evolution, the main difference be- 
tween the MOND gravity and the Newtonian gravity with 
dark matter is the self-gravity of the disc. Even if the accel- 
eration in MOND scales as M 1 / 2 instead of M in Newtonian 
gravity (BM99), the dark matter halo in the DM model sta- 
bilizes more efficiently the disc. From a given initial state, 
the MOND disc is more unstable than the DM disc in the 
sense that it develops a bar instability sooner, for the same 
Toomre parameter value. 

One of the main effects of the dark matter halo is the dy- 
namical friction experienced by the stellar bar against the 
DM particles. The bar pattern speed is slowed down in the 



DM model. This does not exist in MOND. The bar pattern 
speed in MOND keeps constant all along the evolution, thus 
higher than in the DM model. This has consequences on the 
position of the resonances like corotation. Bar lengths are 
often compared to the corotation radius. In this case, bars 
obtained with MOND end closer to the corotation radius. 

The 3D simulations reveal several differences between 
MOND and Newtonian gravity with dark matter. Peanuts 
are formed in the DM model as well as in the MOND 
model, but peanut lobe positions, which correspond to the 
z-inner Lindblad resonance, depend on the bar pattern 
speed. In MOND, the peanut always remains the same size 
(fif, — est), contrary to the DM model where the lobes 
are radially shifted far from the center (about 12 kpc). In 
MOND, successive instabilities due to self-gravity make the 
vertical velocity dispersion higher, in the outer region of 
galaxies, than in the DM model. There is a higher tendency 
for MOND discs to warp and flare. 

Two mechanisms to weaken a bar have been described. 
First, for galaxies with a massive bulge (early type) a 
peanut resonance can be formed. This vertical motion of 
stars dilutes the bar concentration in the plane and makes 
the bar strength decrease. Secondly, if the disc is cold and 
unstable, it forms a bar so quickly that the orbital structure 
of the bar varies on a time scale shorter than the orbital 
period, and the stars cannot settle on orbit supporting the 
bar. 

The present simulations reveal that the dark matter 
halo has two contradictory influences on the disc stabil- 
ity. On the one hand, the DM halo stabilizes the disc and 
delays the bar formation; on the other hand, it can rein- 
force the bar growth when the bar is forming by accepting 
the angular momentum from the disc stars, in particular 
after the peanut's formation. In contrast, peanut galaxies 
in MOND should have low bar strength. 

Statistically, the MOND bar frequency corresponds bet- 
ter to the observations than to the DM model. Indeed, 
there is a hole in the barred galaxy distribution for low bar 
strength and more galaxies distributed at high bar strength. 
But in this work, only stellar discs are considered with- 
out any gas component. Bar formation and destruction is 
affected by the gas component in the spiral galaxies. In 
particular gas accretion allows galaxies to have several bar 
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cycles (Bournaud & Combes, 2002). Gas components will 
be added in future works. 

Through this work, we help to develop numerical tools 
for testing MOND. Using this code, many physical situ- 
ations could be simulated. More complex systems will be 
studied, such as interacting galaxies where MOND might 
reveal larger differences compared with the DM model. 
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